require(KRLS)
require(ggplot2)
require(grDevices)

# Non-monotonicity
data <- read.csv("data.csv") 

data$protest <- 1*(data$s18f1_lag > 0)
Model <- mirt ~ growthc_t1 + lambda + protest + gdp.gle_t1 + agedem_t1 + nelda19 + legelec + country + factor(year)
Model <- lm(Model, data = data)
y <- model.frame(Model)[,1]
X <- model.matrix(Model)[, -1]

kerreg <- krls(X = X, y = y, whichkernel = "gaussian", lambda = NULL,
               sigma = NULL, derivative = TRUE, binary= TRUE, vcov=TRUE, 
               print.level = 1, L=NULL, U=NULL, tol=NULL, eigtrunc=NULL)

pfun <- function(z, xlab = NULL, ylab = NULL){
plot(X[,z], kerreg$derivatives[,z], col = "grey", pch = 16, ylab = "", xlab = xlab, cex.lab = 1.5, las = 1)
mtext(side=2, ylab, las = 1, cex = 1.5, adj=1.75)
x <- sort(X[,z])
yhat <- predict(loess(kerreg$derivatives[, z] ~ X[, z]), newdata = x, se = TRUE)
lines(x, yhat$fit, type = "l", lwd = 5)
lines(x, yhat$fit - 2*yhat$se.fit, type = "l", lwd = 5, lty = 3)
lines(x, yhat$fit + 2*yhat$se.fit, type = "l", lwd = 5, lty = 3)
abline(h = 0, lty = 2, col="dark red")
}

pdf("../graphs/nonlinear.pdf", width = 10, height = 5)
par(mfrow=c(1, 2), mar = c(5,7,0.1,0.1), mgp = c(2,3/4,.1))
pfun(1, xlab = "x = Lagged economic growth", ylab = expression(paste(frac(partialdiff, paste(partialdiff, x)), hat(y), "(x)")))
pfun(2, xlab = "x = Coup threats", ylab = expression(paste(frac(partialdiff, paste(partialdiff, x)), hat(y), "(x)")))
dev.off()
